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Analysis  of  the  Intrinsic  Mode  Functions  1 

by 

Robert  C.  Sharpley  and  Vesselin  Vatchev 


Abstract 

The  Empirical  Mode  Decomposition  is  a  process  for  signals  which  produces  Intrinsic 
Mode  Functions  from  which  instantaneous  frequencies  may  be  extracted  by  simple 
application  of  the  Hilbert  transform.  The  beauty  of  this  method  to  generate  redundant 
representations  is  in  its  simplicity  and  its  effectiveness. 

Our  study  has  two  objectives:  first,  to  provide  an  alternate  characterization  of  the 
Intrinsic  Mode  components  into  which  the  signal  is  decomposed  and  second,  to  better 
understand  the  resulting  polar  representations,  specifically  the  ones  which  are  produced 
by  the  Hilbert  transform  of  these  intrinsic  modes. 


1  Introduction 

The  Empirical  Mode  Decomposition  (EMD)  is  an  iterative  process  which  decomposes  real 
signals  f{t)  into  simpler  signals  (modes) 


M 

=  (Ll) 

3= 1 

Each  “monocomponent”  signal  ipj  (see  [C95])  should  be  representable  in  the  form 

■0F)  =  r(t)  sin  9(t)  (1.2) 

where  the  amplitude  and  phase  are  both  physically  and  mathematically  meaningful.  Once 
a  suitable  polar  parameterization  is  determined,  it  is  possible  to  analyze  the  function  /  by 
processing  these  individual  components.  Important  information  for  analysis,  such  as  the 
instantaneous  frequency  and  instantaneous  bandwidth  of  the  components,  are  derived  from 
the  particular  representation  used  in  (1.2).  The  most  common  procedure  to  determine  a 
polar  representation  is  the  Hilbert  transform  and  this  procedure  will  be  an  important  part 
of  our  discussion. 

1This  work  was  supported  in  part  by  ONR  Grants  N00014-00-1-0470,  N00014-03-1-0051,  and  NSF  Grant 
DMS-0079549. 

Mathematics  Subject  Classijication  2000.  Primary  94A12,  92C55,  41A58;  Secondary  34B24,  44A15. 
Keywords  and  Phrases.  Intrinsic  mode  function,  empirical  mode  decomposition,  signal  processing,  instan¬ 
taneous  frequency,  redundant  representations,  multiresolution  analysis. 
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In  this  paper  we  study  the  monocomponent  signals  ip,  called  Intrinsic  Mode  Functions  or 
IMF’s,  which  are  produced  by  the  Empirical  Mode  Decomposition,  and  their  possible  repre¬ 
sentations  (1.2)  as  real  parts  of  analytic  signals.  Our  study  of  IMF’s  utilizes  mathematical 
analysis  to  characterize  requirements  in  terms  of  solutions  to  self-adjoint  second  order  or¬ 
dinary  differential  equations.  In  principle,  this  seems  quite  natural  since  signal  analysis  is 
often  used  to  study  complex  vibrational  problems  and  the  processes  which  generate  and 
superimpose  the  signal  components.  Once  this  characterization  is  established,  we  then  focus 
on  the  polar  representations  of  IMF’s  which  are  typically  built  using  the  Hilbert  transform, 
or  more  commonly  referred  to  as  the  analytic  method  of  signal  processing. 

The  difficulty  in  constructing  representations  (1.1)  is  that  the  expansion  must  be  selected 
as  a  linear  superposition  from  a  redundant  class  of  signals.  Indeed,  there  are  infinitely  many 
nontrivial  ways  to  construct  representations  of  the  type  (1.1)  even  in  the  case  that  the 
initial  signal  /  is  itself  a  single  “monocomponent”.  Hence  ambiguity  of  representation,  i.e. 
redundancy,  enters  on  at  least  two  levels:  the  first  in  determining  a  suitable  decomposition 
as  a  superposition  of  signals,  and  the  second,  after  settling  on  a  fixed  decomposition,  in 
appropriately  determining  the  amplitude  and  phase  of  each  component. 

At  the  second  stage,  it  is  common  practice  to  represent  the  component  signal  in  complex 
form 

T(f)  =  r(f)  exp  iQ(t)  (1.3) 

and,  after  a  phase  shift  by  n/2,  to  consider  ip(t)  as  the  real  part  of  the  complex  signal  T,  as 
in  (1.2).  Obviously,  the  choice  of  amplitude-phase  representation  (r,  0)  in  (1.3)  is  essentially 
equivalent  to  the  choice  of  an  imaginary  part  <p : 

r(t)  =  J  ip2  +  (p2  9  it)  =  arc  tan  — ,  (1.4) 

once  some  care  is  taken  to  handle  the  branch  cut.  An  analyzing  procedure  should  produce 
for  each  signal  ip,  a  properly  chosen  companion  0  for  the  imaginary  part,  which  is  unambigu¬ 
ously  defined  and  properly  encodes  information  about  the  component  signal,  in  this  case  the 
IMF.  From  the  class  of  all  redundant  representations  of  a  signal,  once  a  fixed,  acceptable  rep¬ 
resentation,  with  amplitude  r(f)  and  phase  9ft),  is  determined,  the  instantaneous  frequency 
of  'pit)  with  respect  to  this  representation  is  the  derivative  of  the  phase,  i.e.  Oft).  In  this 
case,  a  reasonable  definition  for  the  instantaneous  bandwidth  is  (see  [C95]  for  additional 
motivation).  The  collection  of  instantaneous  phases  present  at  a  given  instant  (i.e.  t  =  f0) 
for  a  signal  /(f)  is  heavily  dependent  upon  both  the  decomposition  (1.1)  and  the  selection 
of  representations  (1.2)  for  each  monocomponent.  The  full  EMD  procedure  is  obviously  a 
highly  nonlinear  process,  which  effectively  builds  and  analyzes  inherent  components  which 
are  adapted  to  the  scale  and  location  of  the  signal’s  features. 

Historically,  there  have  been  two  methods  used  to  define  the  imaginary  part  of  suitable 
signals,  the  analytic  and  quadrature  methods.  The  analytic  signal  method  results  in  a 
complex  signal  that  has  its  spectrum  identical  (modulo  a  constant  factor  of  2)  to  that  of 
the  real  signal  for  positive  frequencies  and  zero  for  the  negative  frequencies.  This  can  be 
achieved  in  a  unique  manner  by  setting  the  imaginary  part  to  be  the  Hilbert  transform 
of  the  real  signal  /(f).  The  Empirical  Mode  Decomposition  (EMD)  of  Huang  ef  al  [H99] 
is  a  highly  successful  method  used  to  generate  a  decomposition  of  the  form  (1.1)  where 
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the  individual  components  contain  significant  information.  These  components  were  named 
Intrinsic  Mode  Functions  (IMF’s)  in  [H99]  since  the  analytical  signal  method  applied  to  each 
such  component  normally  provides  desirable  information  inherent  in  that  mode.  Analytic 
signals  and  the  Hilbert  transform  are  powerful  tools  and  are  well  understood  in  Fourier 
analysis  and  signal  processing,  but  in  certain  common  circumstances  the  analytic  signal 
method  leads  to  undesirable  and  “paradoxical”  results  in  applications  which  are  detailed 
in  [C95].  Results  of  this  paper  provide  further  light  on  the  consequences  of  using  analytic 
signals,  as  currently  applied,  for  estimating  phase  and  amplitude  of  a  signal.  More  results 
along  these  lines  appear  in  [V], 

Section  2  contains  a  brief  description  of  the  EMD  method  and  motivates  the  concept  of 
Intrinsic  Mode  Functions  (IMFs)  which  are  the  main  focus  of  this  paper.  Preliminary  results 
on  self-adjoint  equations  are  reviewed  for  background  for  the  results  that  follow  in  Section  3. 
Section  3  contains  one  of  the  main  results  of  the  paper,  namely,  the  characterization  of 
IMFs  as  solutions  to  certain  self-adjoint  ordinary  differential  equations.  The  proof  involves 
a  construction  of  envelopes  which  do  not  rely  on  the  Hilbert  transform.  These  envelopes 
are  used  directly  to  compute  the  coefficients  of  the  differential  equations.  The  differential 
equations  are  natural  models  for  linear  vibrational  problems  and  should  provide  further 
insight  into  both  the  EMD  procedure  and  its  IMF  components.  Indeed,  signals  can  be 
decomposed  using  the  EMD  procedure  and  the  resulting  IMF’s  used  to  identify  systems  of 
differential  equations  naturally  associated  with  the  components  (see  [V]  for  details). 

The  purpose  of  Section  4  of  the  paper  is  to  further  explore  the  effectiveness  of  the  Hilbert 
analysis  which  are  applied  to  IMF’s  and  to  better  understand  some  of  the  anomalies  that 
are  observed  in  practice.  Examples  are  constructed,  both  analytically  and  numerically,  in 
order  to  illustrate  that  the  assumption  that  an  IMF  should  be  the  real  part  of  an  analytic 
signal  leads  to  undesirable  results.  Well-behaved  functions  are  presented,  for  which  the 
instantaneous  frequency  computed  using  the  Hilbert  transform  changes  sign,  i.e.  the  phase 
is  non-monotone  and  physically  unrealistic.  In  order  to  clarify  the  notions  and  procedures, 
we  briefly  describe  both  analytical  and  computational  notions  of  the  Hilbert  transform. 

Finally  we  take  this  opportunity  to  thank  Professor  Ronald  DeVore  (South  Carolina)  and 
Professor  John  Pierce  (USNA)  for  drawing  our  attention  to  the  Empirical  Mode  Decompo¬ 
sition  and  providing  the  primary  references  [H98,  H99]  from  which  to  proceed. 

2  The  Empirical  Mode  Decomposition  Method 

The  use  of  the  Hilbert  Transform  for  decomposing  a  function  into  meaningful  amplitude 
and  phase  requires  some  additional  conditions  on  the  function.  Unfortunately,  no  clear 
description  of  definition  of  a  signal  has  been  given  to  judge  precisely  whether  or  not  a  function 
is  a  “monocomponent”.  To  compensate  for  this  lack  of  precision,  the  concept  of  “narrow 
band”  has  been  adopted  as  restriction  on  the  data  in  order  that  the  instantaneous  frequency 
be  well  defined  and  make  physical  sense.  The  instantaneous  frequency  can  be  considered 
as  an  average  of  all  the  frequencies  that  exist  at  a  given  moment,  while  the  instantaneous 
bandwidth  can  be  considered  as  the  deviation  from  that  average.  The  most  common  example 
is  considered  to  be  a  signal  with  constant  amplitude,  that  is  r(t )  in  (1.4)  is  a  constant.  Since 
the  phase  is  modulated,  these  are  usually  referred  to  as  frequency  modulated,  or  FM  signals. 
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If  no  additional  conditions  are  imposed  on  a  given  signal,  the  previously  defined  notions 
could  still  produce  “paradoxes”.  To  minimize  these  physically  incompatible  artifacts,  Huang, 
et  al  [H99]  have  developed  a  method,  which  they  termed  the  “Hilbert  view” ,  in  order  to  study 
nonstationary  and  nonlinear  data  in  nonlinear  mechanics.  The  main  tools  which  are  used 
are  the  Empirical  Mode  Decomposition  Method  (EMD)  to  decompose  signals  into  Intrinsic 
Mode  Functions  (IMF’s),  which  are  then  processed  by  the  Hilbert  Transform  to  produce 
corresponding  analytic  signals  for  each  of  the  inherent  modes  . 

In  general,  EMD  may  be  applied  either  to  sampled  data  or  to  functions  of  real  variable 
/(f),  by  first  identifying  the  appropriate  time  scales  that  will  reveal  the  physical  character¬ 
istics  of  the  studied  system,  decompose  the  function  into  modes  if  intrinsic  to  the  function 
at  the  determined  scales,  and  then  apply  the  Hilbert  transform  to  each  of  the  intrinsic 
components. 

In  the  words  of  Huang  and  collaborators,  the  EMD  method  was  motivated  “from  the 
simple  assumption  that  any  data  consists  of  different  simple  intrinsic  mode  oscillations”. 
Three  methods  of  estimating  the  time  scales  of  /  at  which  these  oscillations  occur  have  been 
proposed: 

•  the  time  between  successive  zero-crossings  ; 

•  the  time  between  successive  extrema; 

•  the  time  between  successive  curvature  extrema. 

The  use  of  a  particular  method  depends  on  the  application.  Following  the  development  in 
[H99],  we  define  a  particular  class  of  signals  with  special  properties  that  make  them  well 
suited  for  analysis. 

Definition  2.1  A  function  if(t)  is  defined  to  be  an  Intrinsic  Mode  Function,  or  more 
briefly  an  IMF,  of  a  real  variable  t,  if  it  satisfies  two  characteristic  properties: 

(a.)  if  has  exactly  one  zero  between  any  two  consecutive  local  extrema. 

(b.)  if  has  zero  “local  mean”. 

A  function  which  is  required  to  only  satisfy  condition  (a)  will  be  called  a  weak-IMF. 

In  general,  the  term  local  mean  in  condition  (b)  may  be  purposefully  ambiguous,  but  in 
the  EMD  procedure  it  is  typically  the  pointwise  average  of  the  “upper  envelope”  (determined 
by  the  local  maxima)  and  the  “lower  envelope”  (determined  by  the  local  minima)  of  if. 

The  EMD  procedure  of  [H98]  decomposes  a  function  (assumed  to  be  known  for  all  values 
of  time  under  consideration)  into  a  function-tailored,  hne-to-coarse  multiresolution  of  IMF’s. 
This  procedure  is  extremely  attractive,  both  for  its  effectiveness  in  a  wide  range  of  appli¬ 
cations  and  for  its  simplicity  of  implementation.  In  the  latter  respect,  one  first  determines 
all  local  extrema  (strict  changes  in  monotonicity)  and,  for  an  upper  envelope,  fits  a  cubic 
spline  through  the  local  maxima.  Similarly,  a  cubic  spline  is  fit  through  the  local  minima 
for  a  lower  envelope  and  the  local  mean  is  the  average  of  these  two  envelopes.  (It  is  well 
understood  that  these  are  envelopes  in  a  loose  sense).  If  the  local  mean  is  not  zero,  then 
the  current  local  mean  is  subtracted  leaving  a  current  candidate  for  an  IMF.  This  process 
is  continued  (accumulating  the  local  means)  until  the  local  mean  vanishes  or  is  “sufficiently 
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small"’.  This  process  (inner  iteration)  results  in  the  IMF  for  the  current  scale.  The  ac¬ 
cumulated  local  means  from  this  inner  iteration  is  the  version  of  the  function  scaled-up  to 
the  next  coarsest  scale.  The  process  is  repeated  (outer  iteration)  until  the  residual  is  either 
“sufficiently  small"’  or  monotone. 

In  view  of  the  possible  deficiency  of  the  upper  and  lower  envelopes  to  bound  the  iterates 
and  in  order  to  speed  convergence  in  the  inner  loop,  Huang  et  al  suggest  [H99]  that  the  stop¬ 
ping  criterion  on  the  inner  loop  be  changed  from  the  condition  that  the  “resulting  function 
to  be  an  IMF”  to  the  single  condition  that  “the  number  of  extrema  equals  to  zero-crossings” 
along  with  visual  assessment  of  the  iterates.  This  is  the  motivation  for  our  definition  of 
weak-IMF.  Ideally  in  performing  the  EMD  procedure,  all  stopping  and  convergence  criteria 
will  be  met  and  f(t)  is  then  represented  as 


N 

f(t)  =  ^2^n  +  rN+ 1 

n=l 

where  rN+i  =  /n+i  is  the  residual,  or  carrier,  signal. 

A  primary  purpose  of  the  decomposition  [H99]  is  to  distill  from  a  signal  individual  modes 
whose  frequency  (and  possibly  bandwidth)  can  extracted  and  studied  by  the  methods  from 
the  theory  of  analytic  signals.  More  specifically,  quoting  from  [H99], 

“Having  obtained  the  IMF  components,  one  will  have  no  difficulty  in  applying 
the  Hilbert  transform  to  each  of  these  components.  Then  the  original  data  can 
be  expressed  as  the  real  part (3?)  in  the  following  form: 

f(t)  =  ^^A„,(t)exp^i  J  un{t)  dt 

The  residue  r/v  is  left  on  purpose,  for  it  is  either  a  monotonic  function  or  a 
constant.” 

The  notation  above  uses  un  —  to  refer  to  the  instantaneous  frequency,  where  the  phase 
of  the  n-th  IMF  is  computed  by  9n  :  =  arcta xi(H'ipn/'ipn)- 

2.1  Initial  Observations 

The  first  step  in  a  multiresolution  decomposition  is  to  choose  a  time  scale  which  is  inherent 
in  the  function  f(t)  and  has  relevant  physical  meaning.  The  scales  proposed  in  [H99]  are 
sets  of  significant  points  for  the  given  function  f(t).  Other  possibilities  that  could  be  used 
are  the  set  of  inflection  points  (also  mentioned  by  the  authors),  the  set  of  zero  crossings  of 
the  function  f(t)  —  cos  kt,  k- integer,  or  some  other  characteristic  points. 

The  second  step  is  to  extract  some  special  (with  respect  to  the  already  chosen  time  scale) 
functions,  which  in  the  original  EMD  method  are  called  Intrinsic  Mode  Functions(IMF’s). 
The  definition  of  an  IMF,  although  somewhat  vague,  has  two  parts  (a)  the  number  of  the 
extrema  equals  the  number  of  the  zeros  and  (b)  the  upper  and  lower  envelopes  should 
have  the  same  absolute  value.  As  it  is  pointed  out  in  [H99]  if  we  drop  (b)  we  will  have  a 
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reasonable  (from  practical  point  of  view)  definition  but,  in  the  next  stage,  this  will  introduce 
unrecoverable  mathematical  ambiguity  in  determining  the  modulus  and  phase. 

Therefore  any  modification  of  the  definition  of  IMF  must  include  condition  (a).  The 
practical  implementation  of  the  EMD  uses  cubic  splines  as  upper  ( U(t ))  and  lower  ( L(t )) 
envelopes.  The  nodes  of  these  two  splines  interlace  and  do  not  have  points  in  common.  The 
absolute  value  of  two  cubic  splines  can  be  equal  if  and  only  if  they  are  the  same  quadratic 
polynomial  on  the  whole  data  span,  i.e.  if  the  modulus  of  the  IMF  is  of  the  form  at2+bt+c.  To 
overcome  this  restriction,  we  can  either  modify  the  construction  of  the  envelopes  or,  instead 
of  requiring  U{t)  =  — L(t )  for  all  t,  we  can  require  \U(t)  +  L(t)\  <  e,  for  some  prescribed  e  >  0. 

Recall  that  we  say  a  continuous  function  is  a  weak-IMF  if  it  is  only  required  to  satisfy 
condition  (a)  in  Definition  2.1  of  an  IMF.  One  of  the  main  purposes  of  this  paper  is  to  provide 
a  complete  characterization  of  the  weak-IMF’s  in  terms  of  solutions  to  self-adjoint  ordinary 
differential  equations.  In  a  sense  this  is  natural,  since  one  of  the  uses  of  the  EMD  procedure 
is  to  study  solutions  to  differential  equations  and  vibration  analysis  was  a  major  motivation 
in  the  development  of  the  Sturm-Liouville  theory.  In  the  next  section,  we  list  some  relevant 
properties  of  the  solutions  of  a  self-adjoint  ODE’s  which  will  be  useful  for  our  analysis. 

2.2  Self-adjoint  ODE  and  Sturm-Liouville  systems. 

An  ODE  is  called  self-adjoint  if  can  be  written  in  the  form 

|(pw|)+W)/  =  o.  (2.1) 

for  t  e  (a,  b)  (a  and  b  finite  or  infinite),  where  P  >  0  and  Q  is  continuous.  More  generally 
we  can  consider  a  Sturm-Liouville  equation  (A  real): 

jt  +  (V(t)  -  q(t))f  =  0.  (2.2) 

These  equations  arose  from  vibration  problems  associated  with  model  mechanical  systems 
and  the  corresponding  wave  motion  was  resolved  into  simple  harmonic  waves  (see  [BR]). 


Properties  of  the  solutions  of  self-adjoint  and  Sturm-Liouville  equations 

I.  Interlacing  zeros  and  extrema  If  Q  >  0,  then  any  solution  of  (2.1)  has  exactly  one  maxi¬ 
mum  or  minimum  between  successive  zeros. 


II.  The  Priifcr  substitution  A  powerful  method  for  solving  the  ODE  (2.1)  utilizes  a  transfor¬ 
mation  of  the  solution  into  amplitude  and  phase.  If  the  substitution  P  f  :=  r(t)  cos  6{t)  and 
f(t)  :=  r(t)  sin  6{t)  is  made,  then  the  equation  (2.1)  is  equivalent  to  the  following  nonlinear 
first  order  system  of  ODE’s 


ft  =  «*)  sin2  9  +  ^  cos2  9 

dr  1/1  A 

Tt  -  2  (pw -«(*)) 


(2.3) 

(2.4) 
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Notice  that  if  Q(t)  is  positive,  then  the  first  equation  shows  that  the  instantaneous  frequency 
of  the  IMF’s  is  always  positive,  and  therefore  the  solutions  have  nondecreasing  phase.  The 
second  equation  relates  the  instantaneous  bandwidth  with  P(t),  Q(t )  and  9(t).  The  partial 
decoupling  in  this  form  of  the  equations  is  useful  in  studying  the  behavior  of  the  phase  and 
amplitude. 

III.  The  Liouvillc  substitution  An  ODE  of  the  form  (2.2)  can  be  transformed  to  an  ODE  of 
the  type 

f"  +  (A  —  q(t))f  =  0. 

Moreover,  if  fn(t)  is  a  sequence  of  normalized  eigenfunctions,  then 


fnit) 


nn(t  —  a) 
b  —  a 


+ 


Q(i) 

n 


Additional  properties  of  these  solutions  (e.g.  see  [BR])  suggest  that  the  description  of  IMF’s 
as  solutions  to  self-adjoint  ODE’s  will  lead  to  further  insight. 


3  IMF’s  and  Solutions  of  Self-adjoint  ODE 

In  this  section  we  characterize  weak-IMF’s  which  arise  in  the  Empirical  Mode  Decomposition 
algorithm  as  solutions  of  self-adjoint  ODE’s.  The  main  result  may  be  stated  as  follows. 

Theorem  3.1  Let  f  be  a  real-valued  function  in  C2[a,b],  the  set  of  twice  continuously  dif¬ 
ferentiable  functions  on  the  interval  [a,  b] .  If  both  f  and  its  derivative  f  have  only  simple 
zeros,  then  the  following  three  conditions  are  equivalent: 

(i)  The  number  of  the  zeros  and  the  number  of  the  extrema  of  f  on  [a,  b]  differ  by  at  most 
one; 

(ii)  There  exist  positive  continuously  differentiable  functions  P  and  Q  such  that  f  is  a  solu¬ 
tion  of  the  self-adjoint  ODE 


(P(t)f’(t))’  +  Q(t)f(t)  =  0;  (3.1) 

(iii)  There  exists  an  associated  C2[a,b]  function  h  such  that  the  coupled  system 

f(t)  =  h(t)  =  P(t)f(t),  (3.2) 

holds  for  some  positive  continuously  differentiable  functions  P  and  Q. 

Proof:  We  first  prove  that  condition  (i)  is  equivalent  to  (ii).  That  condition  (ii)  implies  (i) 
follows  immediately  since  Q  is  a  positive  function  and  Property  I  of  the  previous  section 
holds  for  solutions  of  self-adjoint  ODE’s  (see  [BR]). 

The  proof  in  the  opposite  direction  ((i)  implies  (ii))  requires  a  preliminary  result  (see 
Lemma  3.1  that  follows)  on  interpolating  piecewise  polynomials  to  be  used  for  envelopes. 
Let  us  assume  then  that  there  is  exactly  one  zero  between  any  two  extrema  of  /.  For 
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simplicity  we  assume  that  the  number  of  zeros  and  extrema  of  /  on  [a,  b]  are  both  equal  to 
M.  Consider  the  collection  of  ordered  pairs 

(3.3) 

which  will  serve  as  our  knot  sequence.  The  points  satisfy  the  required  interlacing 

condition  (ti  <  Z\  <  t2  <  z2  <  •  •  •)>  where  tj  are  the  extremal  points  for  /  and  Zj  are  its 
zeros.  The  data  a  =  {aj}  are  any  positive  numbers  which  satisfy 

max{|/(fj)|,  |/(tj+i)|}  +  V  <  aj  (3.4) 

for  all  j  =  1,  •  •  • ,  M,  where  77  >  0  is  fixed.  The  following  lemma  provides  a  continuous 
piecewise  polynomial  envelope  for  /  by  Hermite  interpolation. 

Lemma  3.1  Let  f  satisfy  the  conditions  of  Theorem  3.1  and  the  {a.,}  satisfy  the  condi¬ 
tion  (3.4),  then  there  is  a  continuous,  piecewise  quintic  polynomial  R  interpolating  this  data 
with  the  following  properties,  for  all  j: 

(a)  The  extrema  of  R  occur  precisely  at  the  points  tj,Zj; 

(b)  | /|  <  R  with  equality  occuring  exactly  at  the  points  tj. 

(c)  R  is  strictly  increasing  on  ( tj,Zj )  and  strictly  decreasing  on  (zj,tj+ 1). 

(d>  R"(tj)  +  (-iy+7"(c) . 

Proof:  Indeed,  let  the  collection  {aj}  satisfy  (3.4),  where  77  >  0  is  fixed.  Interpolate 
the  data  specified  by  (3.3)  by  a  piecewise  quintic  polynomial  R ,  requiring  in  addition  that 
R'itj)  =  R'(zj)  =  0.  On  each  subinterval  determined  by  the  points  {tj,  Zj},  this  imposes  four 
conditions  on  the  six  coefficients  of  the  local  quintic,  leaving  2  degrees  of  freedom  for  each  of 
the  polynomial  ’pieces’.  Representing  such  a  polynomial  in  its  Taylor  expansion  about  the 
left  hand  endpoint  of  its  interval,  it  is  easy  to  verify  that  we  can  force  that  condition  (c)  hold 
at  each  of  the  knots,  and  that  we  can  require  R"(tj)  >  0.  In  particular,  R  has  its  minima 
at  the  maxima  of  |/|  (i.e.  the  tj)  and  its  maxima  at  the  zeros  of  /  (the  Zj ).  Therefore, 
R"{tj)  >  0  >  (— 1  )j+1f"(tj),  which  verifies  condition  (d)  .  □ 

Remark  3.1  In  general,  any  piece-wise  function  R  constructed  from  functions  < pj(t )  that 
satisfy  the  conditions 


v(yi)  =  vu  <p(y2)  =  v2 

v'ivi)  =  ¥>'(2/2)  =  0,  \<p'(t)\  >  0  for  t  e  (2/1, 2/2) 

will  suffice  in  our  construction.  In  particular,  the  Meyer’s  scaling  function  can  be  used  to 
produce  an  envelope  R  which  satisfies  properties  (a)  and  (b)  of  Lemma  3.1  and  can  be  used 
as  a  basis  for  a  quadrature  calculation  of  instantaneous  phase  (see  [V]).  This  idea  is  implicit 
in  the  development  that  follows. 


Having  constructed  an  envelope  R  for  /,  we  define  the  phase- related  function  S  by 


s(t)  :=  ^  (3.5) 

By  Lemma  3.1,  clearly  |5'(i)|  <  1  for  t  G  [a,  6]  and  \S(t)\  =  1  if  and  only  if  t  —  tj  for 
some  j  =  1,2, ,  M.  Since  /  has  exactly  one  zero  between  each  pair  of  consecutive  interior 
extrema,  then  /,  and  hence  S,  has  alternating  signs  at  the  tj.  Without  loss  of  generality,  we 
may  assume  6{ti)  =  |,  i.e.  t\  is  an  interior  local  maximum,  otherwise  we  could  consider  the 
function  — /  instead  of  /.  Endpoint  extrema  are  easily  handled  separately.  As  we  observed 
during  the  proof  of  Lemma  3.1,  the  function  R  was  constructed  to  be  strictly  increasing 
on  ( tj,Zj )  and  strictly  decreasing  on  (z3,  t3+\).  On  intervals  (tj,tj+ 1),  when  j  is  odd,  the 
function  /  decreases,  is  positive  on  ( tj,Zj ),  and  negative  on  (zj,tj+ 1).  These  properties 
imply  that  S  decreases  on  (tj,  t3+\),  is  positive  on  (t3,z3),  and  negative  on  (z3,  t3+\).  Similar 
reasoning  shows  that  for  j  even,  S  increases  on  (tj,  t3+1),  is  negative  on  (t3,z3),  and  positive 
on  (Zj,tj+ 1). 

Therefore  we  can  represent  S  as 


S(t )  =:  sin  6(t)  (3.6) 

for  an  implicit  function  6(t)  which  satishes  0(tj)  =  ^^-7r  and  9(zj )  =  j  it.  From  these  facts, 
one  easily  checks  that  6  is  a  strictly  increasing  function.  In  fact,  6{t)  will  be  continuously 
differentiable  with  strictly  positive  derivative  on  [a,b\.  To  see  this,  first  recall  that  the 
function  R  has  a  continuous  first  derivative  on  [a,  b],  so  S  is  also  differentiable  and  satishes 

SV)  =  (3.7) 

Therefore  S'  is  continuous  and  by  an  application  of  the  implicit  function  theorem  applied  on 
each  of  the  intervals  (tj,tj+ 1),  6{t)  will  be  continuously  differentiable  with  positive  derivative 
on  each  of  these  intervals.  We  will  apply  L’Hospital’s  rule  in  order  to  verify  the  corresponding 
statement  at  the  extrema  tj.  Differentiate  formally  the  relation  (3.6)  and  square  the  result 
to  obtain  on  each  interval  (t3,  t3+i)  the  identity 

(  m  V  m* 

{)  \cos  {6(t)J  1  ~S2(t)' 


So,  if  T(t)  denotes  the  right-hand  side  of  the  above  relation,  that  is 


T(t)  := 


S'(t)2 
1  -  S2(t)  ’ 


(3.8) 


then  T(t)  is  clearly  continuous  except  at  the  tj  where  it  is  undefined.  We  show,  however, 
that  T  has  removable  singularities  at  these  points.  Both  the  numerator  and  denominator 
are  C2  functions  and  vanish  at  tj,  so  an  application  of  L’Hospital’s  rule  shows 


lim  T(t) 

t—>tj 


lim 

t^tj 


2  S'(t)S"(t) 
-2  S'(t)S{t) 


S"{tj) 
S(tj)  ■ 
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On  the  other  hand,  from  (3.7),  S"(tj )  =  ,  and  so  property  (d)  of  the  previous 

lemma  guarantees  that  this  last  expression  is  strictly  positive.  Hence,  9'  is  a  continuous, 
strictly  positive  function  on  the  interval  [a,b\. 

If  we  use  relations  (3.5)  and  (3.6)  to  write  /  as  /  =  R sin#,  then  a  natural  companion  is 
the  function  h  defined  by 

h{t)  :=  —  R(t)  cos9(t).  (3.9) 

It  follows  from  properties  of  R  and  9  that  h  is  strictly  decreasing  on  (zj,  z3+\  )  when  j  is  odd, 
is  strictly  increasing  on  this  interval  when  j  is  even,  and  has  its  simple  zeros  at  the  points 
tj.  Differentiation  of  (3.9),  provides  the  identity 

h!{t)  =  —R\t)  cos 9{t)  +  R{t)9'{t)  sind(t).  (3.10) 


which  will  be  used  to  complete  the  proof  that  condition  (ii)  of  Theorem  3.1  is  satisfied. 
Indeed,  define  the  functions  P,  Q  appearing  in  equation  (3.1)  by 


Pit)  :  = 


h(t) 

fity 


Qit)  := 


m 

nty 


(3.11) 


From  the  properties  of  h  and  /,  we  see  that  these  are  well  defined,  strictly  positive,  and  with 
continuous  first  derivatives  on  [a,  b ],  except  possibly  at  the  set  of  points  {tj}  and  {zj}.  That 
these  properties  persist  at  these  points  as  well,  we  can  again  apply  L’Hospital’s  rule  and  use 
the  identity  (3.10)  together  with  the  fact  that  9'  is  positive. 

Obviously,  the  equations  (3.11)  are  equivalent  to 


Pit)  f\t)  =  r- h[t ),  Qit)  fit)  =  lift)  (3.12) 


which  in  turn  are  equivalent  to  equations  (3.2).  This  establishes  condition  (ii)  of  Theorem  3.1 
and  also  shows  that  this  condition  is  equivalent  to  condition  (iii).  □ 


Remark  3.2  Observe  that  the  function  h  in  condition  (iii)  of  Theorem  3.1  satisfies  a  related 
self-adjoint  ODE, 

(i)  (P(t)h'(t))'  +  Q(t)h{t)  =  0,  where  P  :=  1/Q  and  Q  :=  1/P  and  P,Q  are  the  coeffi¬ 
cients  of  Theorem  3.1  . 

Moreover,  the  coefficients  P,  Q  satisfy  the  following  conditions: 

(ii)  P,  Q  may  be  represented  directly  in  terms  of  the  amplitude  R{t)  and  phase  9{t)  by 

—  =  9'  +  —  t&n9,  Q  =  9' — —  cotd.  (3.13) 

P  it  it 

(iii)  P,Q  satisfy  the  inequality 

with  equality  iff  R'(t)  =  0  on  [a,b\,  or,  equivalently,  if  f  is  an  FM  signal. 

The  only  statements  in  this  Remark  that  requires  justification  are  equations  (3.13).  These 
follow  directly  by  using  the  Priifer  substitution  in  equations  (3.13):  using  (2.3)  for  9'  and  (2.4) 
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for  R'/R. 


Theorem  3.1  provides  the  desired  characterization  of  weak-IMF’s,  which  we  summarize 
in  the  following  corollary. 

Corollary  3.1  A  twice  differentiable  function  if  on  [a,  b]  is  a  weak-IMF  if  arid  only  if  it 
is  a  solution  of  the  self-adjoint  ODE  of  the  type 

{Pif')'  +  Qif  =  0, 


for  positive  coefficients  P(t),  Q(t). 

If  we  adopt  the  definition  of  an  IMF  given  in  Definitiion  2.1,  then  we  have  a  character¬ 
ization  embodied  in  the  following  statements  summarizing  the  results  and  observations  of 
this  section. 


Theorem  3.2  A  function  if  is  an  IMF  if  and  only  if  it  is  a  weak-IMF  whose  spline  envelopes 
satisfy  the  condition  that  the  absolute  value  of  the  lower  spline  envelope  is  equal  to  the  upper 
envelope  and  this  common  envelope  is  a  quadratic  polynomial.  Furthermore,  the  common 
spline  evelope  is  constant  ( i.e if  is  an  FM  signal)  if  and  only  if  Q(t)  =  1  /Pit)  for  the 
associated  self-adjoint  differential  equation  (3.1). 

The  results  of  this  section  indicate  that  we  can  find  a  meaningful  mathematical  and 
physical  description  of  any  weak-IMF  in  terms  of  solutions  of  self-adjoint  problems.  On 
the  other  hand,  considering  these  as  the  real  parts  of  analytic  signals,  we  show  in  the  next 
section  that  there  exist  functions  if  that  are  IMF’s  satisfying  both  conditions  (a)  and  (b)  of 
Definition  2.1,  but  the  phase  produced  by  using  the  Hilbert  transform  is  not  monotonic,  i.e. 
the  instantaneous  phase  changes  sign. 

4  Example  IMFs  and  the  Hilbert  Transform 

In  this  section  we  analyze  several  examples  that  indicate  the  limitations  of  the  analytic 
method  (i.e.  Hilbert  transform)  to  produce  physically  realistic  instantaneous  frequencies  in 
the  context  of  IMF  analysis.  The  examples  presented  show  that  even  for  some  of  the  most 
reasonable  definitions  for  IMFs  the  Hilbert  Transform  method  will  result  in  instantaneous 
frequencies  which  change  signs  on  intervals  of  positive  measure.  By  a  reasonable  IMF  we 
mean  that  they  satisfy  all  existing  definitions,  including  the  IMF  of  Huang  et  al,  narrowband 
mono-components,  and  visual  tests.  Although  our  examples  are  presented  in  order  to  identify 
possible  pitfalls  in  automatic  use  of  the  Hilbert  transform,  in  the  final  analysis,  practitioners 
in  signal  processing  will  make  the  decision  on  when  the  use  of  analyticity  is  appropriate,  and 
to  what  extent  non-monotone  phase  is  necessary.  We  mention  that  the  examples,  in  some 
sense,  also  provide  a  better  understanding  of  many  of  the  paradoxes  concerning  instantaneous 
phase  and  bandwidth  which  are  detailed  in  Cohen  [C95] . 
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4.1  Hilbert  Transforms 

In  order  to  clarify  the  discussion,  we  begin  with  a  brief  description  of  Hilbert  transforms 
and  analyticity.  In  using  the  terminology  “Hilbert  Transform  method” ,  we  mean  one  of  the 
following 

•  the  conjugate  operator  (or  periodic  Hilbert  transform ): 

i.e.,  the  transform  which  is  defined  for  functions  if  on  the  circle  as  the  imaginary  parts 
of  analytic  functions  whose  real  part  coincides  with  fj,  see  [K,  Z]  for  details.  This 
may  be  identified  with  modifying  the  phase  of  each  Fourier  frequency  component  by  a 
quarter  cycle  delay,  i.e.  the  sgn  Fourier  coefficient  multiplier. 

•  the  continuous  Hilbert  Transform : 

i.e.,  the  transform  for  functions  if  defined  on  the  real  line  which  is  defined  as  the 
restriction  to  of  the  imaginary  part  of  analytic  functions  in  the  upper  half  plane 
whose  real  part  on  3?  is  if.  This  is  well  defined  and  understood,  for  example,  on 
Lebesgue,  Sobolev,  Hardy,  and  Besov  spaces  (1  <  p  <  oo  and  in  certain  cases  when 
p  =  oo).  This  transform  may  be  realized  both  as  a  principal  value  singular  integral 
operator  and  as  a  (continuous)  Fourier  multiplier.  For  details  see  [BS,  Z], 

•  the  discrete  Hilbert  transform: 

i.e.,  a  transform  on  discrete  groups  which  is  applied  to  signals  through  a  multiplier 
operator  of  its  discrete  Fourier  transform.  The  operator  is  computed  by  multiplying 
the  FFT  coefficients  of  a  signal  by  sgn  and  then  inverting.  The  multiplier  may  possibly 
invoke  side  conditions  such  as  those  as  implemented  in  the  built-in  version  of  ‘hilbert’ 
in  Matlab  [M],  We  also  note  that  the  m-hle  ‘hilbert. m’  computes  the  discrete  analytic 
signal  itself  and  not  just  the  imaginary  part. 

In  each  of  these  cases  it  is  clear  that  the  imaginary  part  (in  the  case  of  continuous  functions) 
is  uniquely  defined  up  to  an  arbitrary  numerical  constant  C.  In  Fourier  and  harmonic 
analysis  the  choice  is  usually  made  based  on  consideration  of  the  multiplier  operator  as  a 
bounded  isometry  on  L2.  In  some  of  our  examples,  we  will  consider  functions  on  3?  and 
sample  them  in  order  to  apply  the  Discrete  Hilbert  Transform.  For  periodic  functions  and 
appropriate  classes  of  functions  defined  on  3?,  careful  selection  of  the  sampling  resolution 
(e.g.  Shannon  Sampling  Theorem  [P]  in  the  case  of  analyzing  functions  of  exponential 
type)  will  guarantee  that  sampling  the  continuous  Hilbert  transform  of  the  functions  will  be 
equivalent  (at  least  to  machine  precision)  to  application  of  the  discrete  Hilbert  transform  to 
the  sampled  function.  In  other  words,  these  numerical  operations,  when  carefully  applied, 
will  “numerically  comute” .  It  will  be  clear  if  there  is  a  distinction  between  these  transforms 
and,  from  the  context,  which  one  is  intended. 

One  possible  remedy  in  order  to  try  to  avoid  nonphysical  artifacts  of  the  “analytic” 
method  of  computing  the  instantaneous  frequency  is  to  require  additional  constraints  in  the 
definition  of  an  IMF.  One  such  condition  which  immediately  comes  to  mind  would  be  to 
also  require  an  IMF  to  have  at  most  one  inflection  point  between  its  extrema.  We  show  in 
Example  4.2,  however,  that  even  stronger  conditions  are  still  not  sufficient  to  prevent  sign 
changes  of  the  instantaneous  frequency  when  Hilbert  transforms  are  used  to  construct  the 
phase  and  amplitude  for  a  signal,  that  is,  if  one  considers  an  IMF  as  the  real  part  of  an 
analytic  signal.  In  Propositions  4. 1-4.3  we  consider  the  analytical  properties  of  these  exam- 
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pics  and  show  that  they  are  members  of  large  classes  of  signals  that  behave  similarly  when 
processed  by  the  Hilbert  transform,  or  by  the  computational  Hilbert  transform,  not  matter 
how  finely  resolved.  Finally,  we  conclude  this  section  by  describing  a  general  procedure  that 
adds  a  “smooth  perturbation”  to  well  behaved  signals  and  leads  to  undesirable  behavior  in 
estimating  the  instantaneous  phase.  This  indicates  the  need  for  the  possible  consideration 
of  careful  denoising  of  acquired  signals  before  processing  IMF’s  by  the  Hilbert  method. 

Before  proceeding  it  is  useful  to  briefly  discuss  computational  aspects  of  the  Hilbert 
transform  and  therefore  of  the  corresponding  analytic  signal.  There  are  several  versions 
of  the  discrete  Hilbert  transform,  all  using  the  Discrete  Fourier  Transform  (DFT).  In  the 
study  of  monocomponent  signals  which  are  Fourier-based  and  use  least  squares  norms,  the 
choice  of  the  free  parameter  C  is  normally  chosen  so  that  W'tl’Wp  =  ||iT0lk2>  which  mimics  the 
corresponding  property  for  transforms  on  the  line  and  circle.  As  implemented  by  Matlab, 
however,  it  seems  that  for  many  signal  processing  operations  it  is  preferable  to  choose  the  free 
imaginary  constant  so  that  the  constant  (DC)  term  of  the  signal  is  split  between  the  constant 
and  middle  (Nyquist)  terms  of  the  DFT  of  the  Hilbert  transform.  This  appears  natural 
since  the  Nyquist  coefficient  is  aliased  to  the  constant  term,  see  Marple  [M]  for  details.  An 
additional  side  benefit  of  this  choice  of  C  is  that  it  ensures  that  the  discrete  Hilbert  transform 
will  be  orthogonal  to  the  original  signal,  which  emulates  the  corresponding  property  for  the 
Hilbert  transform  for  the  line  and  circle.  We  note  that  the  discretization  process  does  not 
permit  one  to  maintain  all  properties  of  continuous  versions  of  the  transform  and  some  choice 
on  which  properties  are  most  important  must  be  made  based  on  the  application  area. 

One  serious  numerical  artifact  of  the  computational  Hilbert  transform,  which  typically 
arises  when  it  is  applied  to  non-continuous  periodic  functions,  is  a  Gibbs  effect.  Some  care 
must  be  taken  to  insure  continuity  of  the  (implicitly  assumed)  periodic  signal,  otherwise 
severe  oscillations  will  occur  which  often  mask  the  true  behavior  of  the  instantaneous  fre¬ 
quency.  In  the  examples  considered  in  this  section  the  functions  are  continuous,  although 
in  some  cases  (see  Example  4.2)  the  higher  derivatives  are  not.  In  this  case,  however,  the 
oscillations  due  to  this  lack  of  smoothness  are  minor,  of  lower  order  and  do  not  measurably 
affect  the  computations.  Typically,  we  apply  the  computational  Hilbert  transform  after  the 
supports  of  our  functions  are  rescaled  and  translated  to  the  interval  [— 7r,7r). 

Since  it  may  rightly  be  argued  that  other  choices  of  the  free  parameter  C  in  the  discrete 
Hilbert  transform  may  possibly  alleviate  the  problem  of  nonmonotone  phase,  we  focus  for 
the  most  part  on  examples  for  which  any  choice  of  the  imaginary  constant  in  the  analytic 
signal  (and  consequently  in  the  definition  of  the  discrete  Hilbert  transform)  will  result  in 
undesirable  behavior  of  the  instantaneous  frequencies  obtained  by  the  Hilbert  method.  An¬ 
other  concern  in  computational  phase  estimation  is  how  one  numerically  ‘unwraps’  Cartesian 
expressions  to  extract  phases  for  polar  representations.  We  offer  a  technique  to  avoid  am¬ 
biguous  unwrapping  of  inverse  trigonometric  functions  by  instead  computing  the  ‘analytic’ 
instantaneous  frequency  through  the  formula 

=  ip(t)Hip'(t)  -  (Hip(t)  +  C)ip'(t) 

cU'  +  cy  + 

where  9c  is  the  phase  corresponding  to  a  given  choice  of  the  constant  C.  We  use  this  identity 
throughout  to  compute  instantaneous  frequencies  for  explicitly  defined  functions  ip  which  are 
either  periodic  or  defined  on  the  line.  Discrete  versions  using  first  order  differences  are  also 
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suitable  for  computing  instantaneous  phase  for  discretely  sampled  signals.  The  identity  (4.1) 
follows  by  implicitly  differentiating  the  expression  tan  Qq  =  (Hif  +  C)/ if  and  using  the  fact 
that  the  Hilbert  transform  is  translation-invariant. 

We  end  this  subsection  with  an  general  observation  concerning  the  application  of  the 
Hilbert  transform  to  IMF’s,  which  follows  from  Theorem  3.1. 


Corollary  4.1  Suppose  that  if  is  a  periodic,  weak-IMF  and  T  is  the  corresponding  analytic 
function  with  imaginary  part  the  conjugate  operator  Hif.  If  (r,  6)  are  the  corresponding  an¬ 
alytic  amplitude  and  phase  for  the  pair  (if,  Hif),  then  the  coefficients  ( P ,  Q )  of  an  associated 
differential  equation  (2.1)  determined  by  a  Priifer  relationship  (3.13)  must  satisfy 


Hip 

if 


Hif 
if '  ’ 


(4.2) 


whenever  these  two  expressions  make  sense.  In  particular,  a  necessary  and  sufficient  condi¬ 
tion  that  the  coefficients  (P,  Q )  of  the  ODE  be  positive  (i.e.  a  physcially  reasonable  vibrational 
system),  is  that  Hif  should  be  positive  exactly  where  if  increases,  and  if  should  be  positive 
exactly  where  Hif  is  increasing. 


Proof:  This  follows  immediately  from  Theorem  3.1  and  the  Priifer  representation  of  the 
coefficients  which  is  given  in  equation  (3.13).  □ 


4.2  Example  IMFs 

The  first  examples  of  IMF’s  we  wish  to  consider  are  a  family  of  27r-periodic  functions  which 
have  the  property  that  the  conjugate  operator  and  the  discrete  Hilbert  transform  (applied 
to  a  sufficiently  refined  sampling)  differ  only  by  the  addition  of  an  imaginary  constant. 

Example  4.1  Let  e  be  a  real  parameter.  We  consider  the  family  of  continuous  2i r  periodic 
functions 

ife(t)  :=  eecos(t)  sin(e sin(t)).  (4.3) 

Observe  that  the  Hilbert  transform  ofife  is  Hife(t )  =  — eecos^  cos(e  sin(t))  +  C,  where  the 
constant  C  is  a  free  parameter  that  one  may  choose.  In  fact,  the  analytic  signal  T  with  real 
part  ife  is  unique  up  to  a  constant  and  may  be  written  as 

Te(f)  =  —ieeeit  +  iC. 

For  particular  values  of  e  the  function  can  be  used  as  a  model  of  signals  with  interesting 
behavior.  For  example,  ife  for  e  <  2.9716  is  an  FM  function  and  on  any  finite  interval  the 
number  of  the  zeros  differs  from  the  number  of  extrema  by  a  count  of  at  most  one. 

As  one  particular  example  of  the  Hilbert  method  for  computing  instantaneous  phase  for 
IMF’s,  we  fix  in  (4.3)  the  special  choice  of  eo  =  2.97  and  set 

if(t)  :=ifeo(t).  (4.4) 
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Figure  1:  ipe,  an  IMF  with  poor  Hilbert  transform. 


The  graph  of  ip( t )  is  shown  in  Fig  3.1.  In  Proposition  4.1  below,  we  show  that  "0(f)  is  an 
IMF  according  to  the  definition  in  [H98] ,  but  for  any  choice  of  the  constant  C  in  the  Hilbert 
transform,  the  instantaneous  frequency  obtained  from  the  corresponding  analytic  signal  Teo 
changes  sign. 

We  first  verify  the  corresponding  fact  in  the  case  of  discrete  signals.  We  sample  ip( t ) 
uniformly  with  increment  A  =  7t/128  (vector  length  =  1024)  on  the  interval  [— 47t,47t  —  A], 
The  graph  of  the  Hilbert  transform  and  corresponding  instantaneous  frequency  of  ip  obtained 
by  using  Matlab’s  built-in  “hilbert.m”  function  are  shown  in  Figure  2,  parts  (a)  and  (b), 
respectively.  We  mention  that  for  this  data  the  choice  of  constant  chosen  by  Matlab  to 
meet  its  criteria  is  C  =  1.  Although  other  choices  for  C  may  decrease  the  intervals  of 
non-monotonicity  of  the  phase,  the  artifact  will  persist  for  all  choices. 

The  next  proposition  shows  that  the  computational  observation  using  the  discrete  Hilbert 
transform  is  a  consequence  of  the  continuous  transform  and  cannot  be  corrected  by  other 
choices  of  the  imaginary  constant  or  by  a  finer  sampling  rate. 

Proposition  4.1  The  function  ip  defined  by  (f-3)  is  an  IMF  in  the  sense  of  [H98],  but  its 
instantaneous  frequency  computed  by  the  Hilbert  transform  (with  any  choice  of  imaginary 
constant  C)  changes  its  sign  on  any  interval  of  length  at  least  n. 

Proof.  We  first  show  that  ip  is  a  weak-IMF.  Clearly  ip  is  27r-periodic  and  an  odd  function 
and  so  we  only  need  to  consider  it  on  the  interval  [0,  it).  The  first  derivative  of  ip  is 

ip\t)  =  e0eeocos{t)  cos (t  +  e0  sin(t))  (4.5) 

and  is  zero  iff  u(t)  :=  t  +  e0sin(f)  =  =^-tt  for  some  integer  k.  Since  v'(t)  =  1  +  e0  cos(t) 
has  exactly  one  zero  Zq  in  [0, 7r)  ( cos  is  monotone  in  [0, 7r) ) ,  the  function  u(t)  is  increasing 
on  [0,  Zq),  decreasing  on  (£0,  tt) ,  with  end  values  z/(0)  =  0  and  v{i r)  =  it.  To  show  that  ip  has 
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a.  b. 

Figure  2:  the  Hilbert  method  for  ipe:  a.  Imaginary  part  of  discrete  analytic  sig¬ 
nal;  b.  Instantaneous  frequency. 


only  one  extremum  on  [0, 7 r),  it  suffices  to  show  that 

v(zo)  <  (4-6) 

since  then  the  only  extremum  of  ip  on  [0,  n)  will  be  the  point  e«  where  i /(ejw)  =  f .  At  the 
point  z0  we  have  cos(zo)  =  —  l/e0,  which  implies  both  7t/2  <  z0  <  tt  and 

sin(^o)  =  V1  ~  (V eo)2- 
Hence  from  the  dehnition  of  u  it  follows  that 

v{zQ)  =  z0  +  \Jtl~  1- 

This  implies  that  the  condition  (4.6)  is  equivalent  to  z0  <  fvc  —  ^/cq  —  1.  But  cos  is  negative 
and  decreasing  on  [7t/2,  7t],  so  we  see  that  the  desired  relationship  (4.6)  just  means  that 
cos(^o)  >  cos(§7r  —  a/cq  —  1)  should  hold.  The  numerical  value  of  the  expression  on  the  right 
is  smaller  than  —0.3382,  while  cos(zo)  =  —  1/eo  >  —0.3368,  hence  the  condition  (4.6)  holds 
and  "0  has  exactly  one  local  extremum  in  [0, 7 r).  Finally,  since  eo  <  7r,  the  only  zeros  of  ip  are 
clearly  at  the  endpoints  0  and  7r,  which  verifies  that  ip  is  an  weak-IMF. 

To  see  that  ip  is  in  fact  an  IMF,  we  need  to  verify  the  condition  on  the  upper  and  lower 
envelopes.  Recall  that  it  is  2n  periodic  and  odd,  therefore  it  has  exactly  one  minimum  in  the 
interval  [— 7T,0].  The  cubic  spline  fit  of  the  maxima  (upper  envelope)  will  be  the  constant 
function  identically  equal  to  1.  Similarly  the  cubic  spline  interpolant  of  the  minima  (lower 
envelope)  will  have  constant  value  —1.  This  persists  even  for  sufficiently  large  intervals  if 
one  wishes  to  take  finitely  supported  functions.  Hence  the  function  ip  satisfies  the  envelope 
condition  for  an  IMF  from  [H98] .  We  note  that  the  general  proof  to  show  that  for  each  0  < 


16 


e  <  e,  is  an  IMF  follows  in  a  similar  manner,  where  e  is  the  solution  to  the  transcendental 
equation  1/e  =  sin(v/e2  —  1)  which  arises  in  the  limiting  cases  above.  We  observe  that 
e  ~  2.9716. 

Next  we  prove  that  for  any  selection  of  constant  C,  the  corresponding  instantaneous 
frequency  9'c  for  ^  which  is  derived  from  an  analytic  method  through  (4.1)  will  have  nontrivial 
sign  changes.  The  denominator  in  formula  (4.1)  is  always  positive  so  it  will  suffice  to  prove 
that  the  numerator  of  9'c  changes  sign  for  any  choice  of  C.  Using  (4.3),  (4.5)  and  the  formula 

=  e0  ee°  cos ®  sin (t  +  e0  sin (t)) 
we  can  simplify  the  numerator  of  9’c  to  the  expression 

e0  ee°  cos ^  (e£°  cos<^  cos(t)  —  C  cos (t  +  e0  sin(i)))  , 
and  so  the  sign  of  the  term  inside  the  parentheses 

Nt(C)  :=  e£°cos(t)  cos (t)  -  C cos(t  +  e0sin(t)) 

determines  the  sign  of  9'c(t )  at  each  point  t  e  [0, 7r).  First  observe  that  Nt(C)  is  a  linear 
function  of  C  for  fixed  t.  For  each  value  of  C  there  is  a  point  in  [0,  n)  at  which  9'c  is  negative, 
in  fact  Ni  g(C)  <  —0.04  for  C  <  40  while  N0A(C)  <  —8  for  C  >  30.  Similarly  for  any  value 
of  C  there  is  a  point  at  which  9'c  is  positive  since  N0A(C)  >  4  for  C  <  13  and  Nx  (C)  >  2 
for  C  >  0.  By  continuity  we  see  that  for  each  value  of  the  constant  C  the  instantaneous 
frequency  9'c  obtained  via  the  Hilbert  transform  is  positive  and  negative  on  sets  of  positive 
measure.  □ 

Finally,  we  mention  that  the  L2  bandwidth  of  the  analytic  signal  T  corresponding  to  a 
signal  also  depends  on  the  choice  of  the  imaginary  constant  C.  If  T  is  written  in  polar 
coordinates  as  4 /(t)  =  r(t)eld^\  the  average  frequency  (cu)  and  the  bandwidth  u2  have  been 
defined  in  [C95]  as  the  quantities 

(uj)  = 


where  S(u)  is  the  spectrum  (Fourier  transform)  of  ^{t).  The  second  equation  in  the  dis¬ 
played  sequence  (4.8)  follows  immediately  from  Plancherel’s  theorem  along  with  standard 
properties  of  the  Fourier  transform.  If  one  chooses  the  constant  C  in  the  Hilbert  transform 
so  that  HV’lh  =  then  the  computed  bandwidth  is  v2  =  0.1933  with  mean  frequency 
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<  u)  >=  2.7301.  The  discrete  Hilbert  transform  computed  by  matlab  for  the  sampled  ip  has 
the  same  L2  bandwidth  and  mean  frequency. 

Summarizing,  we  observe  that  the  example  ip  given  in  (4.4)  is  a  function  which  is 

(i)  an  IMF  in  the  sense  of  [H98] ; 

(ii)  a  monocomponent  in  the  sense  of  [C95],  i.e.  its  L2  bandwidth  is  small; 

(iii)  “visually  nice” , 

but  the  analytic  method  fails  to  produce  a  monotone  phase  function. 

Remark  4.1  The  example  considered  in  Proposition  4-1  also  shows  that  adding  the  require¬ 
ment  that  the  Hilbert  transform  (with  a  choice  of  the  additive  constant  C  =  3)  of  an  IMF 
must  also  be  a  weak-IMF,  will  not  be  sufficient  to  guarantee  monotone  phase. 

A  possible  natural  refinement  of  the  definition  of  an  IMF  that  would  exclude  these  func¬ 
tions  from  the  class  of  IMFs  would  be  to  require  in  addition  that  the  first  derivative  of  an 
IMF  also  be  an  IMF,  or  at  least  that  the  number  of  the  inflection  points  equals  the  number 
of  extrema  to  within  a  count  of  one  (i.e.  a  weak-IMF).  The  next  example  of  a  damped 
sinusoidal  signal  (i.e.  an  amplitude  modulated  signal )  shows  that  restrictions  along  these 
lines  will  not  be  able  to  avoid  the  same  problem  with  instantaneous  frequencies.  We  note 
that  this  particular  signal  ip  is  considered  in  [H98],  but  for  the  range  of  t  from  1-512  sec. 
Since  the  function  is  not  continuous  periodic  over  this  range,  the  expected  Gibb’s  effect  at 
the  points  t  —  1  and  512  appears  in  that  example,  but  is  absent  here. 

Example  4.2  Let  ip(t )  =  exp(— O.Olf)  cos -^nt,  8  <  t  <  520,  then  ip  is  a  continuous  func¬ 
tion  (of  period  32)  with  a  discontiniuty  in  the  first  derivative  at  t  =  8.  The  signal  ip  and  all 
its  derivatives  are  weak-IMFs.  Both  the  conjugate  operator  and  the  computational  Hilbert 
transform  ( applied  to  the  sampled  function  with  At  =  1)  result  in  a  sign  changing  instanta¬ 
neous  frequency  for  any  choice  of  the  constant  C .  In  Figure  3,  we  have  provided  a  plot  of 
ip  and  the  optimal  instantaneous  frequency  (over  all  possible  C)  which  is  computed  by  the 
Hilbert  transform  method.  The  values  of  both  the  continuous  and  computational  results  are 
to  within  machine  precision  at  the  plotted  vertices. 

In  order  to  verify  the  properties  of  this  example,  we  proceed  as  earlier  in  Example  4.1 
by  first  verifying  the  corresponding  fact  in  the  case  of  discrete  signals.  We  sample  ip(t) 
uniformly  with  increment  A  =  0.1  (vector  length  =  5121)  on  the  interval  [8,  520].  The  graph 
of  the  Hilbert  transform  and  corresponding  instantaneous  frequency  of  ip  obtained  by  using 
Matlab’s  built-in  “hilbert”  function  are  shown  in  Fig  3.3,  parts  (a)  and  (b),  respectively. 
The  next  proposition  shows  that  although  other  choices  of  the  constant  C  may  decrease  the 
interval  where  the  instantaneous  frequency  is  negative  there  is  no  value  for  C  for  which  it  is 
nonnegative  on  [8,520].  Analogous  to  Example  4.1,  it  can  be  shown  that  the  instantaneous 
frequency  changes  its  sign  for  any  choice  of  the  constant  C. 

Proposition  4.2  The  function  ip  in  Example  4.2  and  all  its  derivatives  are  weak-IMFs 
on  the  interval  [8,  512]  whose  instantaneous  frequencies  computed  by  the  Hilbert  transform 
method  change  sign  for  any  choice  of  C. 
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a.  b. 

Figure  3:  Graphs  for  Example  4.2:  (a)  the  IMF;  (b)  its  instantaneous  frequency. 


Proof.  To  simplify  the  notation,  we  denote  by  ip  the  function  in  Example  4.2  and  use  ip 
to  denote  ip  under  the  required  linear  change  of  variable  from  [8,  520]  to  [— tt,  7t]  in  order  to 
apply  the  continuous  Hilbert  transforms  to  the  periodic  function.  In  this  case,  ip  will  be  of 
the  form 

ip(r)  =  c  exp(o:r)  sin  (kr) 

where  k  =  16.  Next  note  that  the  derivatives  of  ip  are  all  of  a  similar  form:  ip^n\t)  = 
c\eat  cos(kt  +  C2).  In  particular,  each  derivative  is  just  a  constant  multiple  of  ip  with  a 
constant  shift  of  phase  and  hence  are  weak-IMF’s  for  any  n  —  0,1,...,  00.  From  formula  (4.1) 
applied  at  the  zeros  of  ip  we  have  0'c(zj )  =  —  where  the  Hilbert  transform  is  dehned 

through  the  conjugate  operator  (see  [Z,  K] )  represented  as  a  principal  value,  singular  integral 
operator 

Hip(zj)  =  —  p.v.  [  ip(t)  cot(— - )  dt.  (4.9) 

^  J  —  7T  2 

The  standard  identity 

t  1 

sm(kt)  cot(-)  =  (1  +  cos(kt ))  +  2  cos  (it)  (4-10) 

£=1 

from  classical  Fourier  analysis  permits  us  to  evaluate  this  expression  by 

HiP(zj)  =  -  [  eatT[t )  dt,  (4.11) 

^  J  —  IT 

where  T  is  an  even  trigonometric  polynomial  of  degree  16  with  coefficients  depending  on 
the  Zj.  Hence  the  values  of  Hip  at  the  zeros  of  ip  can  be  evaluated  exactly  (with  Maple  for 
example).  For  Zi  :=  |7T  and  Z2  '■=  j|vr,  the  corresponding  values  of  the  Hilbert  transform  may 
be  estimated  by  Hip{z\)  <  0.051  and  Hip(z2)  >  0.073  which  verifies  that  Hip(z2)  >  Hi[>(zi). 
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On  the  other  hand,  if'(zi)  =  —keazi  <  0  and  'tp'(z2)  =  keaz 2  >  0,  therefore  0'{zf)  = 
~ H4i- i)+c  negahive  for  C  <  —Hfi{zf)  <  —0.051  and  0'(z2)  =  —  is  negative  for 

C  >  —H'fi^zf)  >  —0.073.  From  the  fact  that  —H'tf^zf)  <  —Hif(zi)  we  conclude  that  for  any 
C  the  instantaneous  frequency  is  negative  for  at  least  one  of  the  points  Z\  or  z2-  Finally,  for 
the  extrema  of  0,  say  t  =  £  we  have  #'(£)  =  cFdif'  ,  where  c  is  a  positive  constant  for 

any  choice  of  C  and  it  is  easy  to  verify  that  there  exists  a  value  £  such  that  F70'(£)0(£)  >  0. 
Hence  #'(£)  >  0  for  any  choice  of  C.  □ 

We  observe  that,  under  the  relaxed  condition  allowing  the  difference  between  the  upper 
and  lower  envelopes  to  be  within  a  given  tolerance,  0  and  its  derivatives  up  to  some  finite 
order  are  (strong)  IMF’s  and  the  computational  Hilbert  transform  method  produces  a  nar¬ 
row  bandwidth  approximately  equal  to  0.0625. 


The  next  result  provides  general  information  about  the  behavior  of  the  instantaneous 
frequency  O'  from  any  polar  representation  of  0(f)  =  r(f)  sin  Oft)  in  terms  of  a  relation 
between  the  amplitude  r  and  0. 


Lemma  4.1  Suppose  that  0  is  a  weak-IMF,  r(t )  >  0  is  an  amplitude  such  that  0(f)  = 
r(t)  cos  Oft).  Further,  suppose  that  at  some  point  t  =  r,  |0(r)|  0  rfr)  and  0(r)  /  0  .A 
necessary  and  sufficient  condition  for  O'  [r )  to  vanish  is  that 


0'(r)  r'(r) 

0(r)  r(r) ’ 

that  is,  that  the  logarithmic  derivative  of  ^  should  vanish  at  t  =  r. 
Proof.  Since  r  >  0  we  can  differentiate  the  relation  cos  6  =  -  and  get 


(4.12) 


— sin(0)  = 


0V  —  00  0/0'  r' 


r  \  0 


(4.13) 


To  prove  necessity,  suppose  that  6'{r)  =  0.  Then  since  0(r)  0  0.  it  follows  from  the  identity 
(4.13),  that  = 

To  prove  sufficiency  it  is  enough  to  notice  that  in  the  event  the  left-hand  side  of  (4.13) 
vanishes  at  f  =  r,  but  0'{r)  f  0,  then  sind(r)  must  vanish.  Hence  |cos6l(r)|  =  1  or 
|0(r)|  =  r(f),  which  is  a  contradiction.  Hence  0'{t)  =  0.  □ 


Looking  back,  one  can  see  that  Lemma  4.1  can  be  used  to  motivate  the  proof  of  the  char¬ 
acterization  theorem  for  weak-IMF’s  (Theorem  3.1).  Indeed,  for  the  envelopes  constructed 
there,  0  and  (-  were  forced  to  have  different  signs  and  therefore  they  can  not  be  equal  at 
any  point  that  is  not  a  zero  of  0.  From  Lemma  4.1,  it  follows  that  O'  does  not  change  sign 
between  any  two  zeros  of  0.  Since  O'  is  continuous  and  was  forced  to  be  nonzero  at  the  zeros 
of  0,  we  have  that  O'  cannot  change  sign. 

Proposition  4.3  Let  va(t)  be  an  even  function  defined  on  (— 7r,7r]  such  that  ua(0)  =  0  and 
eVa/\\eVa  10!  — >  5,  where  5  is  the  Dirac  delta  function.  Define  0a(f)  :=  cos (kt).  Then 
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there  exists  a  value  of  a0  sufficiently  large  such  that  the  analytic  instantaneous  frequency  for 
ip a  for  any  a  >  a0  changes  sign  for  any  choice  of  the  constant  used  in  defining  the  Hilbert 
transform  Hipa. 


Proof:  Recall  from  (4.9)  that  a  Hilbert  transform  of  ipa  at  a  point  t  G  (— 7T,  7 r]  is  Hipa(t)  +  C, 
where  C  is  an  arbitrary  real  constant  and 


HlPaft) 


i  r  ,  ,  N  t-r  , 

p.v.-  /  cot  — - —  dr 


n 


(4.14) 


is  the  conjugate  operator  for  periodic  functions.  Using  two  applications  of  the  identity  (4.1), 
we  observe  that  the  analytic  method  produces  an  instantaneous  frequency  of  the  form 


,,  yjyK  -  m.  +  cx 

c  (H4>a  +  cy-+i>i 


FW'„-CLii, 


(4.15) 


where  R  and  L  are  positive  functions  on  (— 7r,7r].  Let  Zj  =  —k  +  1  <  j  <  k  be  the 

zeros  of  ipa,  then 

sgnfip'fiz.j))  =  {-If  (4.16) 

and  by  (4.15),  with  C  —  0,  it  follows  that 


y'a(zj) 

Hlpa(Zj) 


(4.17) 


and  consequently 

sgn(0'o(zj))  =  (-1  )J+1sgn(H'ipa(zj)).  (4.18) 

The  proof  of  the  lemma  will  be  completed  if  we  can  show  that  for  sufficiently  large  a 
there  is  an  index  J  for  which  two  consecutive  values  of  Htpa(zj)  have  the  same  sign 


sgn  (H'ipa(zj))  =  sgn  (. Hipa{zJ+1 ))  =:  a. 


(4.19) 


When  C  —  0  this  follows  immediately  from  equation  (4.18).  For  C  0,  we  use  the  analogue 
of  (4.18), 

sgn(0'c )  =  - sgn(ip'a )  sgn(Hfia  +  C )  (4.20) 

which  follows  immediately  from  (4.15).  In  the  case  sgn(C )  =  a,  this  last  identity  shows 
that  6'c  has  different  signs  at  the  endpoints  of  (zj,zj+f),  since  ip'a  does.  For  the  final  case, 
sgn(C)  =  —a,  we  observe  that  Hf>a  and  ip'a  are  both  odd  functions  since  ipa  is  even.  By 
considering  —  zj  and  —  zj+\  in  place  of  zj  and  z,j+  \ ,  we  see  that  sgnH{gpa )  =  —a  =  sgn(C)  at 
these  two  points  and  so  once  again  from  (4.20)  ,0c  has  different  signs  at  the  endpoints.  Hence 
by  the  continuity  of  6'C)  there  are  nonempty  intervals  where  the  instantaneous  frequency  takes 
on  opposite  sign. 

Therefore  to  complete  the  proof,  we  must  verify  (4.19),  i.e.,  for  parameter  a  >  0  suffi¬ 
ciently  large,  there  is  a  pair  of  consecutive  points  zj,  zj+ 1,  such  that  Hpipa )  does  not  change 
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sign.  Evaluating  the  conjugate  operator  at  the  zeros  x  =  zj  in  (4.9),  we  can  proceed  as  in 
Proposition  4.2  using  the  periodicity  of  ifa  and  the  trigonometric  identity  (4.10),  to  obtain 


Hlfa(Zj) 


c 


c 


c 


eVa(t)  CQS^f^  CQt  "'■>  fit 

eVa(t+zj)  sin(fci)  cot  L  dt 

eMt+*i)pk(t)  dt, 


where  in  the  last  identity  Pk(t)  is  a  trigonometric  polynomial  of  degree  k.  Therefore  it  follows 
that 


lim 

a— >00 


Hlfa(Zj) 


Pk(  0)  =  -cot^.. 


holds.  For  any  zm,  zm+\  G  (0, 7 r)  it  is  follows  easily  that  there  exists  a  >  0  such  that 
sgn(Hifa(zm ))  =  sgn(Hifa(zm+ 1)).  Hence  for  a  is  sufficiently  large,  ifa  is  a  weak-IMF.  □ 


We  note  that  the  arguments  in  Proposition  4.3  can  also  be  used  to  explain  the  behavior 
of  9'  in  Example  4.2. 

Example  4.3  We  illustrate  in  Figure  4  the  use  of  Proposition  4-3  in  producing  additional 
weak-IMF' s  with  non-monotone  phase.  For  the  sample  function  if,  we  set  k  —  16  and  let  va 
be  a  gaussian  with  standard  deviation  s  =  .01  and  centered  at  the  origin.  The  perturbation 
is  applied  at  both  t\  —  0  and  t2  =  7t/32.  The  function  if  is  displayed  in  part  (a),  its  hilbert 
transform  in  part  (b),  and  its  instantaneous  frequency  in  part  (c). 

For  this  same  function,  in  Figure  4(d)  we  illustrate  the  application  of  Lemma  4-1 ■  The 
instantaneous  frequency  changes  sign  when  the  logarithmic  derivative  of  -  vanishes  at  points 
other  than  at  an  acceptable  zero:  either  a  zero  of  (i)  if  or  of  (ii)  its  Hilbert  transform,  i.e. 
points  where  \  f\  =  r.  Notice  that  the  endpoints  of  the  two  intervals  where  the  instantaneous 
frequency  becomes  negative  corresponds  precisely  to  the  four  (non- acceptable)  zeros  of  the 
logarithmic  derivative  ofif/r 


Example  4.4  An  informative  example  of  function  which  may  be  considered  a  true  IMF  is 
given  by  the  function 


if  it)  =  (f2  T  2)  cos(7rsin(8t))/16,  —  47t  <  t  <  47t  (4-21) 

which,  along  with  its  instantaneous  frequency,  is  plotted  in  Figure  5.  Notice  that  t2  +  2  may 
be  regarded  as  an  envelope  of  if  and  that  it  is  close  but  different  than  the  upper  envelope 
produced  by  a  cubic  spline  fit  through  the  maxima.  Recall  from  the  observation  in  Section  2.1 
that  a  necessary  and  sufficient  for  the  IMF  envelope  condition  (b)  of  Definiton  2.1  to  be 
satisfied  is  that  those  envelopes  reduce  to  a  quadratic  polynomial.  This  example  shows  that 
the  sifting  convergence  criterium  in  the  EMD  process  for  measuring  the  difference  of  the 
absolute  values  of  the  upper  and  lower  envelopes  should  be  chosen  with  care.  If  fact,  it  can  be 
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Sample  IMF  from  Example  4.3:  a.  Cosine  signal  with  strong  pei 
bation  at  0  and  7t/16;  b.  the  Hilbert  tranform  of  ip  near  the  pe: 
bations;  c.  the  instantaneous  frequency  of  ip  near  the  perturbat 
d.  the  logarithmic  derivative  of  ip/r  near  the  perturbations. 


wwwwwww\MAAAzl'liii  hj 
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Figure  5:  Plot  of  the  example  IMF  defined  by  equation  (4.21):  a.  The  IMF  if 
with  a  parabolic  envelope,  b.  Instantaneous  frequency. 


easily  verified  that  smoothing  off  the  endpoint  data  of  this  example  will  result  in  a  function 
for  which  the  EMD  residual  can  be  made  visually  negligible  after  a  single  sifting.  There  are 
many  variations  of  these  examples  to  produce  similar  behavior. 

We  conclude  this  section  with  a  procedure  that  adds  a  smooth  disturbance  at  an  ap¬ 
propriate  scale  to  any  reasonable  function  in  such  a  way  that  the  function  maintains  its 
smoothness  and  analytical  profile,  that  is,  no  additional  extrema  are  introduced  and  the 
existing  extrema  are  only  perturbed,  but  the  resulting  function  has  non-monotone  analytic 
phase.  By  a  reasonable  function  if,  we  will  mean  an  IMF  in  the  strongest  sense,  which  we 
call  a  Hilbert-IMF. 

Definition  4.1  A  function  is  called  a  Hilbert-IMF  if  it  satisfies  the  following  conditions: 
(i.)  if  is  an  IMF  in  the  sense  of  the  definition  in  [H98]; 

(ii.)  the  analytic  signal  of  if  (i.e.  via  the  Hilbert  transform)  T  =  reld ,  r'  and  9'  are 
smooth  functions  and  9'  >  0; 

(iii.)  the  weighted  L2  bandwidth  is  small. 


The  idea  of  the  perturbation  procedure  is  based  on  the  fact  that  by  multiplying  the 
analytic  function  T  =  reie  by  another  analytic  function,  say  T  =  rye*01,  the  result  TT  = 
rriei(0+9i)  js  also  analytic  with  analytic  amplitude  rri  and  analytic  phase  9  +  9\.  Since  9' 
and  9\  are  smooth  functions,  in  order  to  force  the  instantaneous  frequency  of  TT  to  change 
sign  it  suffices  to  choose  T  such  that  9[(T)  <  — 9'{T )  at  some  point  T.  One  way  we  can 
insure  that  the  weighted  L2  bandwidth  of  TT  remains  small  is  to  localize  the  perturbation  T 
to  a  small  interval  I,  i.e.  both  ry  and  9\  should  decay  rapidly  to  zero  outside  /.  Further,  to 
guarantee  that  the  real  part  of  TT  is  an  IMF  in  the  sense  of  [H98],  the  added  perturbation 
must  only  result  in  a  small  deviation  of  the  zeros  and  extrema  of  the  original  IMF  T,  and 
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should  not  introduce  additional  zeros  nor  extrema.  This  is  achieved  by  incorporating  an 
additional  tuning  parameter  a  into  our  perturbation  TCT,  for  small  values  of  a. 

The  technique  used  in  the  proof  of  the  results  below  shows  that  there  are  many  functions 
that  can  be  used  as  such  perturbation  functions.  We  consider  one  particular  smooth  function 
that  is  constructed  from  the  modified  Poisson  kernel 


Vp(t)  ■= 


(1-P)2 

1  —  2 p  cos  t  +  p2  ’ 


0  <  p  <  1. 


(4.22) 


It  is  well  known  that  its  conjugate  function  is  Hyp(t)  =  l-^coTt+y2  •  Although  it  is  an 
abuse  of  notation,  we  will  refer  to  these  simply  as  y  and  Hy  with  the  understanding  that 
the  parameter  p  is  implicitly  present.  We  define  the  perturbation  T  in  terms  of  y  =  yp  by 


T  =  Tp  :=  exp(—Hy  +  iy). 


(4.23) 


and  observe  that  as  the  parameter  p  approaches  1  from  below,  the  function  T p  becomes  very 
localized. 

The  perturbed  IMF  is  set  to  3?(TT<T)  =  re~aHy  008(6*  +  cry)  for  some  real  a.  The  idea  in 
brief  is  to  select  a  small  and  p  sufficiently  close  to  1  so  that  the  change  in  the  functional  values 
are  also  small,  i.e.,  the  zeros,  extrema,  and  the  extremal  values  are  perturbed  slightly  from 
the  original  IMF.  On  the  other  hand,  the  corresponding  instantaneous  frequency  becomes 
9'  +  cry'  (see  Lemma  4.2  below).  Moreover,  y'  has  one  local  minimum  that  is  negative  with 
magnitude  depending  on  p.  In  the  special  case  when  r  =  eAt  and  9  =  mt  on  an  interval  (the 
length  of  the  interval  can  be  arbitrarily  small),  we  prove  in  Corollary  4.2  that  there  exists 
a  subinterval  and  values  of  the  parameters  a  and  p  such  that  under  mild  conditions,  the 
perturbed  function  satisfies  all  the  properties  i-iii),  but  has  non- monotone  phase.  From  the 
proof  and  by  continuity  it  is  then  clear  that  the  new  instantaneous  frequency  can  be  made 
negative  on  an  interval  while  preserving  all  other  properties  listed  in  i-iii).  A  similar  result 
holding  for  more  general  functions  is  established  in  Corollary  4.3. 

To  show  that  the  perturbed  function  is  a  weak-IMF  we  utilize  the  logarithmic  derivative 
as  p  — >  1“  and  the  following  technical  Lemma,  where  we  establish  that  the  maximum  of 
Hy'  and  its  value  at  the  minimum  of  y'  behave  asymptotically  as  a  finite  multiple  of  the 
minimum  value  of  y' . 


Lemma  4.2  Let  0  <  p  <  1  and  y  =  yp  be  defined  as  in  equation  (4-22),  then  the  following 
properties  hold: 


(a) 

(b) 

(c) 


for  all  a  >  0  the  function  T0'  defined  in  equation  (4-23)  is  analytic  with  amplitude 
e~aHy  and  phase  ay; 

if  —p  :=  min y'  =  y'(t0),  then  lim  — —  —  =  —=. 

i-  p  V3 


lim 

p-i- 


P 

max  |  Hy’ 


3^ 

“s- 


(d)  The  function  Hy'  is  even  with  exactly  one  positive  zero,  tz  which  satisfies  0  <  to  <  tz 

and  lim  V-OlL  —  o. 
p^i-  /i 
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Proof.  Part  (a)  follows  from  the  construction  of  the  analytic  function  T  and  the  fact  that 
|T|  >  0.  To  establish  part  (b)  we  determine  the  minimizer  of  y' ,  which  we  denote  f0,  from 
the  equation  y"(to )  =  0,  which  is  equivalent  to  the  equation 

2pcos2f0  +  (1  +  p 2)  cos t0  —  4p  =  0.  (4.24) 

Hence  there  exists  a  unique  solution  t0  which  satisfies  the  relation  cos  to  =  where 

D  :=  a/p4  +  34p2  +  1.  Substituting  this  explicit  formula  for  cos  to  into  the  expression  for 
/(to),  the  minimum  value  of  y'  can  be  written  as 


/(t0)  =  -2(1  -  pf 


a/2(1  +  p2)D  —  (2p4  +  20p2  +  2) 


(3p2  +  3  —  D)2 

Proceeding  similarly  with  the  expression  for  Hy'(t )  given  by 

TT  uf]  =  2p(!  -p)  (1  +  p2)  cost  —  2p 

y[)  1+p  (1  -  2p  cos  t  +  p2)2 

we  hnd,  after  algebraic  rationalization  and  simplification,  that 

Hy'jtp )  _ 2/2 p _ 

P  /(l  +  p2)£  +  p4  +  10p2  +  l' 


(4.25) 


(4.26) 


(4.27) 


Part  (b)  follows  immediately  by  taking  the  limit  as  p  — »  1_. 

Part  (c)  is  established  in  a  similar  manner  by  observing  from  equation  (4.26)  that 
max\Hy'\  =  Hy'( 0)  =  and  so,  using  the  identity  (4.25),  it  follows  that  con¬ 

verges  to  //  as  p  — >  1_. 

Finally,  for  part  (d)  we  determine  from  (4.26)  that  zeros  of  Hy'  are  exactly  the  roots  of 
the  equation  cos  t  =  Substituting  this  expression  for  cos  tz  into  the  left-handside  of 

the  equation  (4.24)  for  cos  to,  we  get  a  negative  value  —/////  for  the  quadratic  expression 
and  so  cos  tz  <  cos  to,  which  is  equivalent  to  t0  <  tz.  Observing  that  sint2  =  /)/,  we  can 
use  this  identity  to  evaluate  y'{tz)  to  obtain  y  ^  0  as  p  — >  1_.  □ 


In  Corollary  4.2  we  prove  in  the  special  case  r  =  eAt ,  H  <  0  and  d  =  mt  that  we  can 
hnd  values  of  a  and  p  such  that  the  procedure  described  above  produces  a  desired  function 
satisfying  the  properties  i-iii)  but  whose  analytic  instantaneous  frequency  changes  sign.  We 
first  prove  a  milder  version  in  the  following  proposition,  and  then  modify  the  parameter  a 
to  establish  the  stronger  version. 


Proposition  4.4  Let  the  notation  be  as  in  the  previous  lemma  (Lemma  Jh2).  In  particular, 
let  to  be  the  point  which  provides  a  global  minimum  for  if ,  tz  be  the  positive  zero  of  Hy' 
and  p  :=  —  /(t0)-  Suppose  further  that  A  <  0.  ///(t)  :=  exp(At)  cos  mt,  then  there  exist  a 
constant  p*  and  a  point  t*  such  that  for  p*  <  p  <  1 

~  777  777 

/(£)  =  exp{At - Hy(t  —  t *))  cos  (mt  H - y(t  —  t*)),  (4.28) 

p  p 
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is  a  weak-IMF,  but  its  analytic  instantaneous  frequency  vanishes  at  t0  +  t* .  Moreover,  the 
difference  between  absolute  values  of  the  upper  and  lower  cubic  spline  envelopes  is  small 
except  at  the  endpoints  in  the  case  that  A  is  negative. 


~  m 

Proof.  From  the  previous  discussion  and  from  the  representation  if(t)  =  9?(T(f)r (t  —  t*)) 
it  is  clear  that  the  analytic  phase  of  if  is  6  =  mt  +  ^y(t  —  t*).  The  definition  of  p  implies 
that  the  expression  m  +  ^ y'(t  —  t*)  is  nonnegative  and  vanishes  only  at  the  point  to  +  t*, 
hence  6  is  strictly  increasing.  Furthermore,  we  show  that  if  p  is  close  to  1,  the  rapid  decay 
of  y(t  —  t*)  away  from  to  +  t*  will  insure  that  the  zeros  of  if  and  if  are  the  same  in  number 
and  are  separated  only  slightly  from  one  another. 

We  may  assume  that  the  perturbation  y  =  yp  is  added  between  a  maximum  of  if  and  the 
zero  To  immediately  following;  the  other  three  situations  can  be  handled  in  the  same  way 
with  an  appropriate  changes  of  the  signs  of  the  corresponding  expressions.  Denote  by  r_  the 
nearest  point  less  than  To  which  satisfies  tanmr_  =  Any  point  from  the  interval  (r_,  To) 

can  be  picked  for  t*.  We  select  t*  :=  n>+f~ ,  5  :=  'n~r~  and  set  A  :=  (t*  —  8,t*  +  8).  By 
construction  it  is  clear  that  functions  y ,  Hy,  y',  and  Hy'  (translated  by  t*)  tend  uniformly  to 
zero  outside  the  interval  A  as  p  approaches  1”.  Hence  if  uniformly  tends  to  if  outside  A. 

Since  if  and  if'  have  only  simple  zeros,  it  follows  that  there  exists  p\  such  that  for  any 
1  >  p  >  p\  the  perturbed  function  if  is  a  weak-IMF;  even  more,  for  each  zero  and  extrema 

of  if  there  corresponds  exactly  one  zero  and  extrema  of  if.  To  prove  this,  we  consider  the 
functions  on  three  disjoint  sets,  a  subinterval  A*  of  A  (to  be  determined),  the  set  A\A*, 
and  the  complement  of  A. 

We  first  consider  the  set  of  values  t  in  the  complement  of  A.  Assume  that  there  is  a 
sequence  of  p's  approaching  1  from  below  so  that  there  is  an  extrema  of  if,  say  re,  such  that 
in  a  neighborhood  of  that  extrema  if  has  at  least  three  extrema  (the  extrema  must  be  odd  in 
number  since  if  tends  uniformly  to  if).  By  Rolle’s  theorem  and  the  uniform  convergence  as 
p  — >  1_,  it  follows  that  if'  has  a  multiple  zero  at  re,  which  is  a  contradiction.  Hence  outside 
A,  if  is  an  IMF. 

Consider  now  the  interval  A.  To  follow  the  changes  of  the  extrema  of  if,  we  consider  the 
logarithmic  derivative  of  if  given  by 

L  :=  i  =  A  +  -(-Hy'(t  -  t*))  -  (m  +  —y'(t  -  t*))  tan  (mt  +  —y(t  -  t*)). 

if  p  p  p 

We  show  L  is  negative  on  A.  It  then  follows  that  if  has  no  additional  extrema  on  this 

interval,  and  so  if  will  be  an  IMF  in  the  sense  of  the  definition  in  [H98],  but  its  analytic 
instantaneous  frequency  m  +  ^y1  has  a  zero. 

To  see  that  L  is  negative  on  A,  observe  that  for  a  fixed  p  (1  >  p  >  pf)  since  A  <  0, 
tan  (mt  +  ^y(t  —  t*))  >  0,  and  (m  +  ^ y'(t  —  t*))  >  0  it  follows  that  L  <  0  on  the  interval 
A*  :=  (t*  —  tz,t*  +  tz),  where  tz  is  specified  in  part  (d)  of  Lemma  4.2.  From  the  proof 
of  that  Lemma,  we  see  that  tz  approaches  0  as  p  approaches  1“  and  hence  we  can  pick  p2 
(1  >  p2  >  pi)  such  that  A*  C  A  for  each  p  for  which  1  >  p  >  p2- 
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On  the  other  hand,  for  any  t  outside  A*,  we  have  that  m  +  y y'(t )  >  m  +  y y'(tz ) 
and  from  Lemma  4.2(d)  it  follows  that  there  exists  1  >  p3  >  p2  such  that  the  inequality 
m  +  — y'(t )  >  y  holds  for  any  1  >  p  >  p3  and  any  t  G  A\A*.  Finally  from  Lemma  4.2, 

we  can  pick  1  >  p*  >  p3  such  that  the  inequality  max  dO  I  <  AA  holds  for  any  p  >  p*.  The 
choice  of  the  point  t*  and  the  fact  that  y  is  a  positive  function  provide  the  inequality 
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ta n(mt  H - y{t  —  t  ))  >  ta n(mr_)  = 

for  any  t  e  A\A*.  Using  the  above  estimates  and  the  assumption  A  <  0  we  have  that 

~  maxlF/VI  m  .  m  .  ...  „ 

L  <  m - tan(mi  H - y(t  —  t  ))  <  0, 

p  2  p 

for  any  t  e  A\A*  and  any  1  >  p  >  p*,  which  completes  the  proof.  □ 

Corollary  4.2  Let  if  6e  the  Hilbert-IMF  considered  in  Proposition  4-4 ■  Then  there  exist 
a  >  3^  such  that  the  small  perturbation  of  ip  =  3?(vhrCT),  satisfies  parts  (i)  and  (in)  of 
the  definition  of  a  Hilbert-IMF,  but  has  does  not  have  monotone  analytic  phase. 

Proof.  Since  and  all  other  related  functions  considered  in  Proposition  4.4  depend  contin¬ 
uously  on  a,  and  for  o  =  —  —  we  have  that  9'  vanishes  only  at  the  point  t0  +t*,  it  follows 
that  any  increase  of  a  forces  the  instantaneous  frequency  to  be  negative  in  a  neighborhood 
of  that  point.  On  the  other  hand,  by  the  choice  of  p*(  1  >  p*  >  0)  from  Proposition  4.4 
and  the  uniform  convergence,  it  follows  that  there  exists  a*  >  —  ^  such  that  the  perturbed 

function  is  a  weak-IMF  outside  A  and  La*  is  negative  on  A;  i.e.  there  are  no  additional 

zeros  and  extrema  on  A.  Hence  if  is  a  nicely  behaved  function  with  analytic  instantaneous 
frequency  that  changes  sign  on  an  interval  of  positive  measure.  □ 


Remark  4.2  The  condition  A  <  0  can  be  relaxed  to  A  <  which  agrees  with  the  estimate 
in  Lemma  4 -2(b).  The  proof  of  the  Proposition  4-4  w'dh  that  restrictiion  requires  further 
technical  estimates  as  in  Lemma  4-%(d)  for  a  point  tv  such  that  0  <  t0  <  tv  and 


lim 

p-i- 


y'ifv) 

p 


rj  —  1 


for  a  fixed  0  <  rj  <  1.  Details  of  the  estimates  are  similar  so  we  do  not  include  them  here. 


Example  4.5  To  illustrate  the  above  construction,  we  consider  the  function  4'  =  e4lt,  —it  < 
t  <  it,  and  apply  the  procedure  twice  for  p  =  0.95  and  a  =  0.31,  once  at  the  point  t\  =  —2.1 
and  again  at  the  point  0.2.  The  resulting  signal  and  its  analytic  instantaneous  frequency  are 
shown  on  Figure  6. 


For  any  nice  function  (e.g.  a  Hilbert-IMF)  if  =  r  cos  9  with  only  simple  zeros,  it  is  clear 
from  the  identity  —  —  —  tan  6  9'  that  if  there  exist  a  zero  for  which  r'  0  then  the 

logarithmic  derivative  y-  and  the  instantaneous  bandwidth  y  have  the  same  sign  in  a  one¬ 
sided  neighborhood  of  that  zero.  Then  if  the  perturbation  is  added  in  that  neighborhood, 
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Figure  6:  Plot  of  the  IMF  in  Example  (4.5):  a.  an  IMF  perturbed  by  a  smooth 
perturbation;  b.  Instantaneous  frequency. 


the  proof  of  Proposition  4.4,  without  significant  modifications,  can  be  used  to  prove  similar 
results  so  long  as  y  +  max/  9'  is  negative  in  a  neighborhood  of  that  zero,  as  established 
in  the  next  Corollary.  We  note  that  by  choosing  p  closer  to  1,  the  interval  /  can  be  made 
arbitrarily  small. 

Corollary  4.3  Let  ip  =  r  cos  6  be  the  restriction  to  the  circle  of  any  function  analytic  in  a 
disk  of  radius  larger  than  one.  Assume  that  ip  is  an  IMF  in  the  sense  of  [H98]  with  amplitude 
r  and  monotone  phase  9  which  are  defined  using  the  Hilbert  transform.  Suppose  further  that 
ip  has  only  simple  zeros,  that  there  exists  a  zero  z0  of  ip  for  which  ^  and  have  the  same 

sign  in  a  one-sided  neighborhood  I  of  z0,  and  that  y  +  max/  9'  <  0  on  I ,  then  there  exist 

parameters  p*  and  a* ,  and  a  point  t*  such  that  the  function  ip(t )  =  K(4'(t)r<T*  [t  —  t *))  is  an 
IMF  with  zeros  and  extrema  which  are  close  perturbations  of  those  of  ip,  but  with  an  analytic 
instantaneous  frequency  9'{t)  —  9'  +  cry' it  —  t*)  which  changes  sign. 

Proof.  As  in  Proposition  4.4,  it  is  enough  to  prove  that  9'{T )  =  0  for  some  point  T .  We 
may  assume  from  the  hypothesis  that  9'  >  0,  then  the  modified  instantaneous  frequency 
9 [t)  =  9'{t)  +  cry' it  —  t*)  >  0  for  small  a.  Hence  if  a  is  continuously  increased,  by  continuity 
we  will  reach  a  value  G\  for  which  9'(T)  =  0  for  some  point  T  and  is  positive  elsewhere.  If  p 
is  then  increased  close  to  1,  and  is  adjusted  accordingly,  we  can  localize  the  perturbation 
on  an  arbitrarily  small  interval  with  9'  vanishing  at  a  point. 

Notice  that  the  instantaneous  bandwidth  does  not  change  sign  on  an  interval  that  contains 
both  /  and  z0  as  an  interior  point.  We  may  assume  that  L  <  0  on  /,  then  the  logarithmic 
derivative  of  ip  is 

~  T1  (t)  ~  ~ 

L(t )  =  — — —  crHy'(t  —  t*)  —  9\t)  tan  9{t). 

r{t) 
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From  the  choice  of  a  it  follows  that  max j  O'  —  ay  >  0  and  hence  a  <  ma*jg  ,  For  the  new 
instantaneous  bandwidth  as  p  — ■>  1_,  we  have  y  —  aHy'(t  —  t*)  <  y  +  max/  6'  "iaxTy  I  — >  rL  y 
-^=  max/  9'  <  0  on  I,  where  the  last  inequality  is  the  assumption  on  if  relating  instantaneous 

bandwidth  and  frequency.  Since  /  is  a  finite  interval,  there  exists  p*  such  that  L  <  0  for  any 
p  satisfying  1  >  p  >  p*.  On  the  other  hand,  since  9'  >  0  and  tan  6  >  0  on  /  it  follows  that 
L  <  0  on  I.  All  other  steps  are  the  same  as  in  the  proof  of  Proposition  4.4.  □ 

Remark  4.3  The  condition  on  the  logarithmic  derivatives  (of  the  function  if  and  its  am¬ 
plitude  r)  to  have  the  same  sign  in  a  neighborhood  of  a  zero  is  equivalent  to  r'(zo)  ^  0,  or 
r'(z0 )  =  0  but  r'lf'if  >  0  on  I.  In  other  words  if  all  other  requirements  are  met,  the  procedure 
works  more  generally  than  in  the  case  of  envelopes  considered  in  the  Theorem  3.1. 

Remark  4.4  The  procedure  for  adding  perturbations  to  a  nice  function  can  be  also  used  for 
removing  certain  types  of  noise.  If  a  function  has  a  negative  instantaneous  frequency  on 
some  small  interval,  then  we  can  apply  the  procedure  with  a  perturbation  r0-  using  negative  a 
in  order  to  remove  negative  instantaneous  frequencies,  but  still  preserve  the  general  features 
(zeros,  local  extrema)  and  smoothness  class  of  the  original  function. 


References 

[BS]  C.  Bennett  and  R.  Sharplcy,  Interpolation  of  Operators,  Academic  Press,  New  York, 

1988. 

[BR]  G.  Birkhoff  and  G.C.  Rota,  Differential  Equations,  (4-th  ed)  Wiley,  New  York,  NY, 

1989. 

[C95]  L.  Cohen,  Time-Frequency  Analysis,  Prentice-Hall  Signal  Processing  Series,  Upper 
Saddle  River,  NJ,  1995. 

[H98]  N.E.  Huang,  Z.  Shen,  S.R,  Long,  M.C.  Wu,  H.H.  Shih,  Q.  Zheng,  N.C.  Yen,  C.C. 
Tung,  and  H.H.  Liu,  The  empirical  mode  decomposition  and  the  Hilbert  spectrum 
for  nonlinear  and  non-stationary  time  series  analysis,  Proc.  Royal  Soc.  Lond.  A  454 
(1998),  903-995. 

[H99]  N.E.  Huang,  Z.  Shen,  and  S.R.  Long,  A  new  view  of  nonlinear  water  waves:  the  Hilbert 
spectrum,  Annu.  Rev.  Fluid  Mech.  31  (1999),  417-457. 

[K]  Y.  Katznelson,  An  Introduction  to  Harmonic  Analysis,  (2-nd  ed.),  Dover  Pubns,  1976. 

[M]  S.L.  Marple,  Computing  the  discrete-time  analytic  signal  via  FFT,  IEEE  Transactions 
on  Signal  Processing  47  (1999),  2600-2603. 

[P]  M.A.  Pinsky,  Introduction  to  Fourier  Analysis  and  Wavelets,  Brooks/Cole,  Pacific 
Grove,  2002. 


30 


[S]  E.M.  Stein,  Singular  Integrals  and  Differentiability  Properties  of  Functions,  Princeton 
University  Press,  Princeton,  1970. 

[V]  V.  Vatchev,  Intrinsic  Mode  Functions  and  the  Hilbert  Transform,  Ph.D  Dissertation, 
Department  of  Mathematics,  University  of  South  Carolina,  2004. 

[Z]  A.  Zygmund,  Trigonometric  Series,  Cambridge  University  Press,  Cambridge,  1969. 


Industrial  Mathematics  Institute 
Department  of  Mathematics 
University  of  South  Carolina 
Columbia,  SC  29208 

E-mail:  sharpley@math.sc.edu,  vatchev@math.sc.edu 


31 


